/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: AliAnalysisTaskESDMuonFilterO.cxx 57113 2012-06-12 08:52:42Z fca $ */

//
// Add the muon tracks to the generic AOD track branch during the 
// filtering of the ESD. 
//
// Authors: R. Arnaldi 5/5/08 and L. Aphecetche January 2011
//
// Note that we :
//   - completely disable all the branches that are not required by (most) the muon analyses,
//     e.g. cascades, v0s, kinks, jets, etc...
//   - filter severely the tracks (keep only muon tracks) and vertices (keep only primary -including
//     pile-up - vertices) branches 
// 
// (see AddFilteredAOD method)
//

#include "AliAnalysisTaskESDMuonFilterO.h"

#include "AliAODDimuon.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliAODExtension.h"
#include "AliAODMCParticle.h"
#include "AliAODMuonReplicator0.h"
#include "AliAODVertex.h"
#include "AliAnalysisFilter.h"
#include "AliAnalysisManager.h"
#include "AliCodeTimer.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDMuonTrack.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMultiplicity.h"
#include <TChain.h>
#include <TFile.h>
#include <TParticle.h>
#include "AliESDtrackCuts.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
#include "TF1.h"
#include "AliPIDResponse.h"
#include "AliAOD3LH.h"

using std::cout;
using std::endl;
ClassImp(AliAnalysisTaskESDMuonFilterO)
ClassImp(AliAnalysisNonMuonTrackCuts)

////////////////////////////////////////////////////////////////////////

AliAnalysisNonMuonTrackCuts::AliAnalysisNonMuonTrackCuts()
{
  // default ctor 
}

Bool_t AliAnalysisNonMuonTrackCuts::IsSelected(TObject* obj)
{
  // Returns true if the object is a muon track
  AliAODTrack* track = dynamic_cast<AliAODTrack*>(obj);

  ULong_t  status;

  if(track){
    
    status  = (ULong_t)track->GetStatus();

    if(track->GetTPCNcls() > 80 &&
       track->Chi2perNDF() < 5  &&
       track->IsOn(AliAODTrack::kTPCrefit) &&
       track->IsOn(AliAODTrack::kTPCin)    &&
       !track->IsOn(AliAODTrack::kITSpureSA))
      {
	return kTRUE;
      }
  } 
  
  else 
    return kFALSE;
  

}

AliAnalysisNonPrimaryVertices::AliAnalysisNonPrimaryVertices()
{
  // default ctor   
}

Bool_t AliAnalysisNonPrimaryVertices::IsSelected(TObject* obj)
{
  // Returns true if the object is a primary vertex
  
  AliAODVertex* vertex = dynamic_cast<AliAODVertex*>(obj);
  if (vertex)
    {
      if ( vertex->GetType() == AliAODVertex::kPrimary     ||
	   vertex->GetType() == AliAODVertex::kMainSPD     ||
	   vertex->GetType() == AliAODVertex::kPileupSPD   ||
	   vertex->GetType() == AliAODVertex::kPileupTracks||
	   vertex->GetType() == AliAODVertex::kMainTPC )
	{
	  return kTRUE;
	}
    }
  
  //  enum AODVtx_t {kUndef=-1, kPrimary, kKink, kV0, kCascade, kMulti, kMainSPD, kPileupSPD, kPileupTracks,kMainTPC};

  return kFALSE;
  
}

AliAnalysisTaskESDMuonFilterO::AliAnalysisTaskESDMuonFilterO(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Int_t nsigmaTrk1,Int_t nsigmaTrk2, Int_t partType1,Int_t partType2):
  AliAnalysisTaskSE(),
  fTrackFilter(0x0),
  fEnableMuonAOD(kTRUE),
  fEnableDimuonAOD(kTRUE),
  fOnlyMuon(onlyMuon),
  fKeepAllEvents(keepAllEvents),
  fMCMode(mcMode),
  fnSigmaTrk1(nsigmaTrk1),
  fnSigmaTrk2(nsigmaTrk2),
  fpartType1(partType1),
  fpartType2(partType2),
  fPIDResponse(0)
{
  // Default constructor
}

AliAnalysisTaskESDMuonFilterO::AliAnalysisTaskESDMuonFilterO(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Int_t nsigmaTrk1,Int_t nsigmaTrk2, Int_t partType1,Int_t partType2):
  AliAnalysisTaskSE(name),
  fTrackFilter(0x0),
  fEnableMuonAOD(kTRUE),
  fEnableDimuonAOD(kTRUE),
  fOnlyMuon(onlyMuon),
  fKeepAllEvents(keepAllEvents),
  fMCMode(mcMode),
  fnSigmaTrk1(nsigmaTrk1),
  fnSigmaTrk2(nsigmaTrk2),
  fpartType1(partType1),
  fpartType2(partType2),
  murep(0),
  fPIDResponse(0)
{
  // Constructor
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::UserCreateOutputObjects()
{
  //-----------------------------------------------
  // Particle Identification Setup (new PID object)
  //-----------------------------------------------
  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  fPIDResponse = inputHandler->GetPIDResponse();
  // Create the output container
  if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::PrintTask(Option_t *option, Int_t indent) const
{
  // Specify how we are configured
  
  AliAnalysisTaskSE::PrintTask(option,indent);
  
  TString spaces(' ',indent+3);
  
  if ( fOnlyMuon ) 
    {
      cout << spaces.Data() << "Keep only muon information " << endl;        
    }
  else 
    {
      cout << spaces.Data() << "Keep all information from standard AOD" << endl;
    }
  
  if ( fKeepAllEvents ) 
    {
      cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;    
    }
  else 
    {
      cout << spaces.Data() << "Keep only events with at least one muon" << endl;
    }
  
  if ( fMCMode > 0 ) 
    {
      cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
    }
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::AddFilteredAOD(const char* aodfilename, const char* title)
{
  // Add an output filtered and replicated aod
  
  AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
  if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");

  AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title);

  if (!ext) return;
  
  if ( fOnlyMuon ) 
    {    
    
      // AliAODMuonReplicator* murep = new AliAODMuonReplicator("MuonReplicator",
      //                                                        "remove non muon tracks and non primary or pileup vertices",
      //                                                        new AliAnalysisNonMuonTrackCuts,
      //                                                        new AliAnalysisNonPrimaryVertices,
      //                                                        fMCMode);
      
      // AliAODMuonReplicator0* murep = new AliAODMuonReplicator0("MuonReplicator0",
      // 							       "remove non interesting tracks",
      // 							       new AliAnalysisNonMuonTrackCuts,
      // 							       new AliAnalysisNonPrimaryVertices,
      // 							       fMCMode);
      if(!murep)delete murep;

      murep = new AliAODMuonReplicator0("MuonReplicator0",
      					"remove non interesting tracks",
      					new AliAnalysisNonMuonTrackCuts,
      					new AliAnalysisNonPrimaryVertices,
      					fMCMode);
      
      
      
      ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
      
      ext->FilterBranch("tracks",murep);    
      ext->FilterBranch("vertices",murep);  
      ext->FilterBranch("dimuons",murep); //per test
      ext->FilterBranch("AliAODVZERO",murep);
      ext->FilterBranch("AliAODTZERO",murep);
      
      cout<<"add filterd aod"<<endl;
      
      if ( fMCMode > 0 ) 
	{
	  // MC branches will be copied (if present), as they are, but only
	  // for events with at least one muon. 
	  // For events w/o muon, mcparticles array will be empty and mcheader will be dummy
	  // (e.g. strlen(GetGeneratorName())==0)
	  
	  ext->FilterBranch("mcparticles",murep);
	  ext->FilterBranch("mcHeader",murep);
	}
    }  
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::Init()
{
  // Initialization
  if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
  //  if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");    
}


//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event					    
  
  Long64_t ientry = Entry();
  if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
  cout<<"--------------->Enter the user exec<------------------------------"<<endl;
  
  //Check the PIDresponse
  // if(!fPIDResponse) {
  //   AliError("Cannot get pid response");
  //   return;
  // }

  //***************************************************

  ConvertESDtoAOD();

  //*************************************************
  
  if(!murep)delete murep;
  
}
//_____________________________________________________________________________
Bool_t AliAnalysisTaskESDMuonFilterO::IsTrackAccepted(AliESDtrack * const track)
{
  
  //! TPC
  
  Int_t    minclsTPC=60;
  //  Int_t    minclsTPC=10;
  Double_t maxchi2perTPCcl=5.;
  
  AliESDtrackCuts *fEsdTrackCuts = new AliESDtrackCuts("esdtrackCuts");
  fEsdTrackCuts->SetRequireITSStandAlone(kFALSE);
  fEsdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
  //  fEsdTrackCuts->SetRequireITSRefit(kTRUE);

  fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
  fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
  fEsdTrackCuts->SetMinNClustersTPC(minclsTPC);
  fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxchi2perTPCcl);

  Bool_t IsTrackAccepted =  fEsdTrackCuts->AcceptTrack(track);
  if (IsTrackAccepted) 
    return kTRUE;
 
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskESDMuonFilterO::IsHighDeDxParticle(AliESDtrack * track) {
  
  //
  // check if particle is Z > 1 
  //

  // if (track->GetTPCNcls() < 60) 
  //   return kFALSE;
  
  Double_t mom = track->GetTPCmomentum();
  //track->GetInnerParam()->GetP(); //solitamente usato
  
  if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
  Float_t dca[2], bCov[3];
  track->GetImpactParameters(dca,bCov);

  //Double_t triggerDeDx = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)7);
  
  
  Double_t triggerDeDx =4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3.4),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
    
    //  if (triggerDeDx > -3 && track->GetTPCsignal() < 1000 /*&& TMath::Abs(dca[0])<3.*/) 
 
  if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<5000 /*&& TMath::Abs(dca[0])<3.*/) 
    return kTRUE;
  
  return kFALSE;
} 

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)
{
  //
  // Setter for the raw PID detector signals
  //

  if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track
    detpid = new AliAODPid();
    SetDetectorRawSignals(detpid,esdtrack);
    aodtrack->SetDetPID(detpid);
  }
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilterO::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)
{
//
//assignment of the detector signals (AliXXXesdPID inspired)
//
  if(!track){
    AliInfo("no ESD track found. .....exiting");
    return;
  }
  
  // TPC momentum
  aodpid->SetTPCmomentum(track->GetTPCmomentum());
  // aodpid->SetTPCTgl(track->GetTPCTgl());

  // aodpid->SetITSsignal(track->GetITSsignal());
  // Double_t itsdedx[4]; // dE/dx samples for individual ITS layers
  // track->GetITSdEdxSamples(itsdedx);
  // aodpid->SetITSdEdxSamples(itsdedx);

  aodpid->SetTPCsignal(track->GetTPCsignal());
  aodpid->SetTPCsignalN(track->GetTPCsignalN());
  
  // if(track->GetTPCdEdxInfo()) 
  //   aodpid->SetTPCdEdxInfo(track->GetTPCdEdxInfo());

  //---fin qui
  
  // //n TRD planes = 6
  // Int_t nslices = track->GetNumberOfTRDslices()*6;
  // TArrayD trdslices(nslices);
  // for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
  //   for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
  // }
 
  // //TRD momentum
  // for(Int_t iPl=0;iPl<6;iPl++){
  //   Double_t trdmom=track->GetTRDmomentum(iPl);
  //   aodpid->SetTRDmomentum(iPl,trdmom);
  // }
  
  // aodpid->SetTRDslices(track->GetNumberOfTRDslices()*6,trdslices.GetArray());
  // aodpid->SetTRDsignal(track->GetTRDsignal());
  
  // //TRD clusters and tracklets
  // aodpid->SetTRDncls(track->GetTRDncls());
  // aodpid->SetTRDntrackletsPID(track->GetTRDntrackletsPID());
  
  // aodpid->SetTRDChi2(track->GetTRDchi2());
  
  // //TOF PID  
  // Double_t times[AliPID::kSPECIES]; track->GetIntegratedTimes(times);
  // aodpid->SetIntegratedTimes(times);
  
  // aodpid->SetTOFsignal(track->GetTOFsignal());
  
  // Double_t tofRes[5];
  // for (Int_t iMass=0; iMass<5; iMass++){
  //   tofRes[iMass]=0; //backward compatibility
  // }
  // aodpid->SetTOFpidResolution(tofRes);
  
}


//______________________________________________________________________________
/*
void AliAnalysisTaskESDMuonFilterO::ConvertESDtoAOD() 
{
  // ESD Muon Filter analysis task executed for each event
  
  AliCodeTimerAuto("",0);
  
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!esd) 
    {
      AliError("Could not get input ESD event");
      return;    
    }
  
  AliMCEventHandler *mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
    
  // Define arrays for muons
  Double_t pos[3];
  Double_t p[3];
  Double_t pid[10];
  AliAODPid* detpid(0x0);
  
  // Primary vertex cut
  
  const AliESDVertex *vtx = esd->GetPrimaryVertexTracks();
  
  if(vtx->GetNContributors()<1) {
    
    // SPD vertex cut
    vtx = esd->GetPrimaryVertexSPD();
    
    if(vtx->GetNContributors()<1) {
      Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
      return; // NO GOOD VERTEX, SKIP EVENT 
    }
  }
  
  if(TMath::Abs(vtx->GetZv())>10) return;
  
  // has to be changed once the muon pid is provided by the ESD
  // for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
  // pid[AliAODTrack::kMuon]=1.;
  
  AliAODHeader* header = AODEvent()->GetHeader();
  AliAODTrack *aodTrack = 0x0;
  AliESDtrack *esdMuTrack = 0x0;
  
  // Access to the AOD container of tracks
  TClonesArray &tracks = *(AODEvent()->GetTracks());
  Int_t jTracks = header->GetRefMultiplicity();
  
  // Read primary vertex from AOD event 
  AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
  if (fDebug && primary) primary->Print();
  
  // Loop on muon tracks to fill the AOD track branch
  Int_t nMuTracks = esd->GetNumberOfTracks();//esd->GetNumberOfMuonTracks();
  cout<<"Numero di tracce: "<<nMuTracks<<endl;
  for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) 
    esd->GetTrack(iTrack)->SetESDEvent(esd);
  

  // // Update number of positive and negative tracks from AOD event (M.G.)
  // Int_t nPosTracks = header->GetRefMultiplicityPos();
  // Int_t nNegTracks = header->GetRefMultiplicityNeg();
  // 
  Int_t nPosTracks = 0;
  Int_t nNegTracks = 0;
  
  // Access to the AOD container of dimuons
  //  TClonesArray &dimuons = *(AODEvent()->GetDimuons());
  
  Int_t nMuons=0;
  Int_t nDimuons=0;
  Int_t jDimuons=0;
  Int_t nMuonTrack[100];
  UChar_t itsClusMap(0);
  Int_t contHe=0;
  Int_t contHelii=0;
  //for(int imuon=0;imuon<100000;imuon++) nMuonTrack[imuon]=0;

  for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
    {
      esdMuTrack = esd->GetTrack(nMuTrack);
      
      //    if (!esdMuTrack->ContainTrackerData()) continue;
      if(!esdMuTrack) continue;
      if(esdMuTrack->Charge()==0) continue;
      if(IsTrackAccepted(esdMuTrack)==kFALSE) continue;
      
      Double_t mom =  esdMuTrack->GetTPCmomentum();
      Float_t dca[2], bCov[3];
      esdMuTrack->GetImpactParameters(dca,bCov);

      if(IsHighDeDxParticle(esdMuTrack)){
	++nMuons;
	contHe++;
	
	cout<<"\n\nPRESO\n\n "<<endl;

	// cout<<"nsigma - 6: "<<fPIDResponse->NumberOfSigmasTPC(esdMuTrack,(AliPID::EParticleType)6)<<endl;
	//	cout<<"nsigma - 7: "<<fPIDResponse->NumberOfSigmasTPC(esdMuTrack,(AliPID::EParticleType)7)<<endl;
	// cout<<"nsigma - 8: "<<fPIDResponse->NumberOfSigmasTPC(esdMuTrack,(AliPID::EParticleType)8)<<endl;
	
      }
      
      UInt_t selectInfo(0);
    
      p[0] = esdMuTrack->Px(); 
      p[1] = esdMuTrack->Py(); 
      p[2] = esdMuTrack->Pz();
      
      pos[0] = esdMuTrack->GetX(); 
      pos[1] = esdMuTrack->GetY(); 
      pos[2] = esdMuTrack->GetZ();
      
      itsClusMap = esdMuTrack->GetITSClusterMap();
      esdMuTrack->GetESDpid(pid);
      
      if (mcH) mcH->SelectParticle(esdMuTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
      
      aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
						    esdMuTrack->GetLabel(), // label
						    p, // momentum
						    kTRUE, // cartesian coordinate system
						    pos, // position
						    kFALSE, // isDCA
						    0x0, // covariance matrix
						    esdMuTrack->Charge(), // charge
						    itsClusMap, // ITSClusterMap
						    pid, // pid
						    primary, // primary vertex
						    kTRUE, // used for vertex fit?
						    //                                              kTRUE,
						    vtx->UsesTrack(esdMuTrack->GetID()), // used for PV?
						    AliAODTrack::kPrimary,// FIXME track type
						    selectInfo); 
      
      Double_t Chi2perNDF=0;
      Int_t  nClustersTPC = esdMuTrack->GetTPCNcls();
    
      if ( nClustersTPC > 5) {
	Chi2perNDF=esdMuTrack->GetTPCchi2()/Float_t(nClustersTPC - 5);
      } else {
	Chi2perNDF=-1.;
      }
      
      aodTrack->SetTPCFitMap(esdMuTrack->GetTPCFitMap());
      aodTrack->SetTPCClusterMap(esdMuTrack->GetTPCClusterMap());
      aodTrack->SetTPCSharedMap (esdMuTrack->GetTPCSharedMap());
      aodTrack->SetChi2perNDF(Chi2perNDF);
      aodTrack->SetTPCPointsF(esdMuTrack->GetTPCNclsF());
      aodTrack->SetTPCNCrossedRows(UShort_t(esdMuTrack->GetTPCCrossedRows()));
      aodTrack->SetFlags(esdMuTrack->GetStatus());
      aodTrack->ConvertAliPIDtoAODPID();
      SetAODPID(esdMuTrack,aodTrack,detpid);//new
      primary->AddDaughter(aodTrack);
    
      if (esdMuTrack->Charge() > 0) nPosTracks++;
      else nNegTracks++;
    
    }

  cout<<"BBB Cont HE3"<<  contHe <<endl;
  
  if(contHe>=1)
    cout<<"SBUNNA"<<endl;

  contHelii=contHelii+contHe;
  
  header->SetRefMultiplicity(jTracks); 
  header->SetRefMultiplicityPos(nPosTracks);
  header->SetRefMultiplicityNeg(nNegTracks);
  header->SetNumberOfMuons(nMuons);
  header->SetNumberOfDimuons(nMuons);
  
  AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
  
  //if ( handler && fEnableMuonAOD && (  (contHe>0)|| fKeepAllEvents ) ) 
  
  if ( handler && fEnableMuonAOD && (contHe>0)) 
    {
     
      cout<<"\n\n\n------------->QUIIIIIIIIIIIIII: "<<contHe<<endl;
      AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
      
      if ( extMuons ){ 
	extMuons->SelectEvent();
	extMuons->Print();
      }
    }
    

  // if ( handler && fEnableMuonAOD && (  (contHe>0) )) 
  // // if ( handler && fEnableDimuonAOD && ( (contHe>1) || fKeepAllEvents )  )
  //   {
  //     cout<<"voglio scrivere i di muoni"<<endl;
  //     AliAODExtension *extDimuons = handler->GetFilteredAOD("AliAOD.Dimuons.root");
  //     if ( extDimuons ) extDimuons->SelectEvent();
  //     cout<<"Ho scritto Di muoni"<<endl;
  //   }

  cout<<"Helii: "<<  contHelii<<endl;
  
}
*/


//---- Funziona


void AliAnalysisTaskESDMuonFilterO::ConvertESDtoAOD() 

{
  // // ESD Muon Filter analysis task executed for each event
  
  // AliCodeTimerAuto("",0);
  
  // AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  // if (!esd) 
  //   {
  //     AliError("Could not get input ESD event");
  //     return;    
  //   }
  
  // AliMCEventHandler *mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
    
  // // Define arrays for muons
  // Double_t pos[3];
  // Double_t p[3];
  // Double_t pid[10];
  // AliAODPid* detpid(0x0);
  
  // // Primary vertex cut
  
  // const AliESDVertex *vtx = esd->GetPrimaryVertexTracks();
  
  // if(vtx->GetNContributors()<1) {
    
  //   // SPD vertex cut
  //   vtx = esd->GetPrimaryVertexSPD();
    
  //   if(vtx->GetNContributors()<1) {
  //     Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
  //     return; // NO GOOD VERTEX, SKIP EVENT 
  //   }
  // }
  
  // if(TMath::Abs(vtx->GetZv())>10) return;
  
  // // has to be changed once the muon pid is provided by the ESD
  // // for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
  // // pid[AliAODTrack::kMuon]=1.;
  
  AliAODHeader* header = AODEvent()->GetHeader();
  AliAODTrack *aodTrack = 0x0;
  AliESDtrack *esdMuTrack = 0x0;
  
  // Access to the AOD container of tracks
  TClonesArray &tracks = *(AODEvent()->GetTracks());
  Int_t jTracks = header->GetRefMultiplicity();
  
  // Read primary vertex from AOD event 
  AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
  if (fDebug && primary) primary->Print();
  
  // Loop on muon tracks to fill the AOD track branch
  //Int_t nMuTracks = esd->GetNumberOfTracks();//esd->GetNumberOfMuonTracks();
  //cout<<"Numero di tracce: "<<nMuTracks<<endl;
  // for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) 
  //   esd->GetTrack(iTrack)->SetESDEvent(esd);
  
  
  // // Update number of positive and negative tracks from AOD event (M.G.)
  // Int_t nPosTracks = header->GetRefMultiplicityPos();
  // Int_t nNegTracks = header->GetRefMultiplicityNeg();
  // 
  Int_t nPosTracks = 0;
  Int_t nNegTracks = 0;
  
  // Access to the AOD container of dimuons
  //  TClonesArray &dimuons = *(AODEvent()->GetDimuons());
  
  Int_t nMuons=0;
  Int_t nDimuons=0;
  Int_t jDimuons=0;
  Int_t nMuonTrack[100];
  UChar_t itsClusMap(0);
  Int_t contHe=0;
  Int_t contHelii=0;
  //for(int imuon=0;imuon<100000;imuon++) nMuonTrack[imuon]=0;

  // for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
  //   {
  //     esdMuTrack = esd->GetTrack(nMuTrack);
      
  //     //    if (!esdMuTrack->ContainTrackerData()) continue;
  //     if(!esdMuTrack) continue;
  //     if(esdMuTrack->Charge()==0) continue;
  //     if(IsTrackAccepted(esdMuTrack)==kFALSE) continue;
      
  //     Double_t mom =  esdMuTrack->GetTPCmomentum();
  //     Float_t dca[2], bCov[3];
  //     esdMuTrack->GetImpactParameters(dca,bCov);

  //     if(IsHighDeDxParticle(esdMuTrack)){
  // 	++nMuons;
  // 	contHe++;
	
  // 	cout<<"\n\nPRESO\n\n "<<endl;

  // 	// cout<<"nsigma - 6: "<<fPIDResponse->NumberOfSigmasTPC(esdMuTrack,(AliPID::EParticleType)6)<<endl;
  // 	//	cout<<"nsigma - 7: "<<fPIDResponse->NumberOfSigmasTPC(esdMuTrack,(AliPID::EParticleType)7)<<endl;
  // 	// cout<<"nsigma - 8: "<<fPIDResponse->NumberOfSigmasTPC(esdMuTrack,(AliPID::EParticleType)8)<<endl;
	
  //     }
      
  //     UInt_t selectInfo(0);
    
  //     p[0] = esdMuTrack->Px(); 
  //     p[1] = esdMuTrack->Py(); 
  //     p[2] = esdMuTrack->Pz();
      
  //     pos[0] = esdMuTrack->GetX(); 
  //     pos[1] = esdMuTrack->GetY(); 
  //     pos[2] = esdMuTrack->GetZ();
      
  //     itsClusMap = esdMuTrack->GetITSClusterMap();
  //     esdMuTrack->GetESDpid(pid);
      
  //     if (mcH) mcH->SelectParticle(esdMuTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
      
  //     aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
  // 						    esdMuTrack->GetLabel(), // label
  // 						    p, // momentum
  // 						    kTRUE, // cartesian coordinate system
  // 						    pos, // position
  // 						    kFALSE, // isDCA
  // 						    0x0, // covariance matrix
  // 						    esdMuTrack->Charge(), // charge
  // 						    itsClusMap, // ITSClusterMap
  // 						    pid, // pid
  // 						    primary, // primary vertex
  // 						    kTRUE, // used for vertex fit?
  // 						    //                                              kTRUE,
  // 						    vtx->UsesTrack(esdMuTrack->GetID()), // used for PV?
  // 						    AliAODTrack::kPrimary,// FIXME track type
  // 						    selectInfo); 
      
  //     Double_t Chi2perNDF=0;
  //     Int_t  nClustersTPC = esdMuTrack->GetTPCNcls();
    
  //     if ( nClustersTPC > 5) {
  // 	Chi2perNDF=esdMuTrack->GetTPCchi2()/Float_t(nClustersTPC - 5);
  //     } else {
  // 	Chi2perNDF=-1.;
  //     }
      
  //     aodTrack->SetTPCFitMap(esdMuTrack->GetTPCFitMap());
  //     aodTrack->SetTPCClusterMap(esdMuTrack->GetTPCClusterMap());
  //     aodTrack->SetTPCSharedMap (esdMuTrack->GetTPCSharedMap());
  //     aodTrack->SetChi2perNDF(Chi2perNDF);
  //     aodTrack->SetTPCPointsF(esdMuTrack->GetTPCNclsF());
  //     aodTrack->SetTPCNCrossedRows(UShort_t(esdMuTrack->GetTPCCrossedRows()));
  //     aodTrack->SetFlags(esdMuTrack->GetStatus());
  //     aodTrack->ConvertAliPIDtoAODPID();
  //     SetAODPID(esdMuTrack,aodTrack,detpid);//new
  //     primary->AddDaughter(aodTrack);
    
  //     if (esdMuTrack->Charge() > 0) nPosTracks++;
  //     else nNegTracks++;
    
  //   }

  // cout<<"BBB Cont HE3"<<  contHe <<endl;
  
  // if(contHe>=1)
  //   cout<<"SBUNNA"<<endl;

  // contHelii=contHelii+contHe;
  
  // header->SetRefMultiplicity(jTracks); 
  // header->SetRefMultiplicityPos(nPosTracks);
  // header->SetRefMultiplicityNeg(nNegTracks);
  // header->SetNumberOfMuons(nMuons);
  // header->SetNumberOfDimuons(nMuons);
  
  AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
  
  //if ( handler && fEnableMuonAOD && (  (contHe>0)|| fKeepAllEvents ) ) 
  
  if ( handler /*&& fEnableMuonAOD && (contHe>0)*/) 
    {
     
      cout<<"\n\n\n------------->QUIIIIIIIIIIIIII: "<<contHe<<endl;
      AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
      
      if ( extMuons ){ 
	extMuons->SelectEvent();
	extMuons->Print();
      }
    }
    

  // if ( handler && fEnableMuonAOD && (  (contHe>0) )) 
  // // if ( handler && fEnableDimuonAOD && ( (contHe>1) || fKeepAllEvents )  )
  //   {
  //     cout<<"voglio scrivere i di muoni"<<endl;
  //     AliAODExtension *extDimuons = handler->GetFilteredAOD("AliAOD.Dimuons.root");
  //     if ( extDimuons ) extDimuons->SelectEvent();
  //     cout<<"Ho scritto Di muoni"<<endl;
  //   }

  //  cout<<"Helii: "<<  contHelii<<endl;

}
//------------------------------------------------
void AliAnalysisTaskESDMuonFilterO::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  //  delete murep;
  
  if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
}
